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Abstract 

The foundations are laid for the numerical computation of the actual worldline for a particle 
orbiting a black hole and emitting gravitational waves. The essential practicalities of this com- 
putation are here illustrated for a scalar particle of infinitesimal size and small but finite scalar 
charge. This particle deviates from a geodesic because it interacts with its own retarded field il) ret . 
A recently introduced [1] Green's function G s precisely determines the singular part, ip s , of the 
retarded field. This part exerts no force on the particle. The remainder of the field = V ,ret — ip s 
is a vacuum solution of the field equation and is entirely responsible for the self-force. A particular, 
locally inertial coordinate system is used to determine an expansion of ip° in the vicinity of the 
particle. For a particle in a circular orbit in the Schwarzschild geometry, the mode-sum decom- 
position of the difference between ip ret and the dominant terms in the expansion of -0 s provide a 
mode-sum decomposition of an approximation for R from which the self-force is obtained. When 
more terms are included in the expansion, the approximation for tp R is increasingly differentiable, 
and the mode-sum for the self-force converges more rapidly. 

PACS numbers: 04.25.-g, 04.20.-q, 04.70.Bw, 04.30.Db 
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I. INTRODUCTION 



In general relativity, a particle of infinitesimal mass will orbit a black hole of large mass 
along a worldline T which is an exact geodesic in the background geometry determined by 
the large mass alone. If the orbiting particle is not infinitesimal, having a small finite mass, 
its orbit will no longer be a geodesic in the background of the larger mass, and gravitational 
waves will be emitted by the system — at infinity. In the neighborhood of the small particle, 
local measurements cannot separately distinguish the background of the large mass from a 
smooth perturbation to it caused by the presence of the smaller mass [1] . The actual orbit of 
the particle can be analyzed in a linearization of the Einstein equations via a perturbation 
expansion in the ratio of the masses. Through first order in this ratio, T is known to be a 
geodesic of a geometry perturbed from the background of the large mass by the presence 
of the smaller one [1]. The difference of the worldline from a geodesic in the background is 
said to arise from the interaction of the orbiting particle with its own gravitational field. It 
is said to result from a "self-force," even though, in the perturbed geometry determined by 
both the small and large masses, the orbit would be observed to be geodesic. 

In a strict sense, if the particle is of infinitesimal size, then its own field is singular 
along its worldline, and there the perturbation analysis fails. This difficulty can be avoided 
by allowing the size of the particle to remain finite while invoking the conservation of the 
stress-energy tensor within a world-tube which surrounds the worldline, in a manner similar 
to Dirac's [2] classical analysis. The balance of energy and momentum indicates how to 
calculate the self-force in a way which is independent of the size of the particle. The limit 
of vanishing size may then be taken without confusion. 

In curved spacetime, analyses beginning with DeWitt and Brehme [3] and subsequently 
by Mino, Sasaki and Tanaka [4] and by Quinn and Wald [5, 6] formally resolve the difficulty 
presented by the singularity in curved spacetime with a Hadamard expansion [3] of the 
Green's function near T. The retarded Green's function G rct (p,p') incorporates physically 
appropriate boundary conditions and describes the field ifj ret of a particle moving through a 
given spacetime. In most past discussions of the self-force, G ret (p,p') is commonly divided 
into a "direct" part, which has support only on the past null cone of the field point p, and 
a "tail" part, which has support inside the past null cone and is a result of the curvature of 
spacetime. The analyses show that ip tml = ip vct — ip dn is necessarily finite at the particle and 
suggest that it is the part of ip which belongs on the right hand side of an equation for the 
self-force [6] 

J = a = qV a ^. (1) 

In these approaches, the use of ?/> tai1 instead of ip ret constitutes a form of regularization of the 
singular ip rct . Actually, while finite, ?/> tai1 is generally not differentiable on the worldline [3] if 
the Ricci scalar of the background is not zero. Similarly, the electromagnetic potential A* ai1 
(respectively, the gravitational metric perturbation h 1 ^ 1 ) is not differentiable at the particle 
if (Rob — ^g a bR)u b (respectively, R ca dbU c u d ) is nonzero in the background. In all such cases, 
some version of averaging must be invoked to make sense of the self-force. Moreover, we 
find it instructive to observe that the tail part of the field is necessarily associated with a 
nonphysical inhomogeneous source, i.e. V a V a, ?/> tai1 7^ 0: cf. Eq. (30). 

In this paper an alternative regularization of the field ip ret is used to compute the self-force 
where, in particular, by regularization we mean not only controlling the singular behavior, 
but also the differentiability. We have recently given a precise procedure for decomposing 
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the retarded field in neighborhood of T into two parts [1] 

where ip s is a solution of the inhomogeneous field equation for the particle, and is determined 
in the neighborhood of the worldline entirely by local analysis via Eq. (36). As both if) ret 
and ip s are inhomogeneous solutions of the same differential equation, it follows that ip R , 
defined by Eq. (2) is necessarily a homogeneous solution and is therefore expected to be 
differentiable on T. In Ref. [1] we showed that ?/> R formally gives the correct self-force when 
substituted on the right hand side of Eq. (1) in place of ^ tai1 . In this paper ?/> R is used for 
an explicit computation of the self-force. We consider ip s to be associated with the Singular 
Source, and ip K with the Regular Remainder. 

While the procedure which follows from Eq. (2) is well understood in principle, its ap- 
plication to physically interesting situations remains a challenge. In this paper we consider 
a particle endowed with a scalar charge q in circular motion about a Schwarzschild black 
hole. On a technical level, a spherical harmonic decomposition of both ip ret and ^> s provides 
the multipole components of each, and the mode by mode sum of the difference of these 
components determines ip R and, thence, the self-force. 

In Section II we give a brief overview of the relation between our work and that of earlier 
authors. We also summarize our analytical results and introduce the additional regularizing 
parameters which allow us to obtain increased convergence in our mode sum representation 
of the self-force. 

A special set of coordinates is described in Section III; the THZ coordinates, introduced 
by Thorne and Hartle [7] and extended by Zhang [8], are locally inertial on a geodesic. 
These coordinates are convenient for describing the scalar wave equation in the vicinity of 
the geodesic, where the metric takes a particularly advantageous form. In Section IV the 
Hadamard expansion for the Green's function, discussed in detail by DeWitt and Brehme 
[3], is described in terms of Synge's [9] "world function" (j(p,p r ), which is defined as half of 
the square of the geodesic distance between two points p and p'. We obtain both cr(p,p') 
and the Hadamard expansion in terms of the THZ coordinates. 

Section V outlines the determination of the regularization parameters given below in 
Eqs. (13) to (15). These results are in agreement with, but extend by going to higher order, 
the work of Barack and Ori [10-12] and Mino, Nakano and Sasaki [12, 13]. 

In Section VI, with a concrete application of our method, we examine a scalar charge in 
a circular orbit of the Schwarzschild geometry at a radius of 10M. It is in this section that 
we see the practical advantage of using a higher order approximation in the regularization 
of ip s . The additional parameters we find enable us to increase dramatically the rate of 
convergence in the self-force summation. 

In several Appendices we include details concerning the THZ coordinates, the mathemati- 
cal analyses which focus on calculation of the regularization parameters and a brief summary 
of details concerning the integration of the scalar wave equation in the Schwarzschild geom- 
etry. 

Notation: a,b,... are four dimensional space-time indices. are three dimen- 

sional spatial indices. (t 8 ,r,0,(f>) are the usual Schwarzschild coordinates. t,x,y,z are 
locally-inertial THZ coordinates attached to the geodesic T, and p 2 = x 2 + y 2 + z 2 . The 
geodesic T is given as x a = z a (r), where r is the proper time along T. The flat spatial 
metric in Cartesian coordinates is The flat Minkowski metric in Minkowski coordinates 
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is r/ab = (— 1, 1, 1, 1). The points p and p' refer to a field point and a source point on the 
world line of the particle, respectively. In the coincidence limit p — > p'. An expression such 
as 0(p n ) means of the order of p n as x l — * in the THZ coordinates. But note that the 
differentiability of such an order term is only necessarily C n_1 at x l = 0. 

II. OVERVIEW OF RELATION TO EARLIER WORK 

Formally, although our approach differs from that of Barack and Ori, our method of 
implementation is similar to that in their pioneering analysis in Ref. [10, 11, 14]. In their 
procedure, which Burko has implemented [15, 16] both for a scalar field with radial and with 
circular orbits of Schwarzschild, the self-force may be thought of as being evaluated from 
[30] 

K CU = lim [J?(p) -J?*®], (3) 

where p' is the event on Y where the self-force is to be determined, p is an event in the 
neighborhood of p', and the relationship between J- a {p) and ip(p) is as given in Eq. (1). To 
make use of this equation, both ^ et (p) and J-^ip) are expanded into multipole £-modes, 
with J-J^ip) determined numerically. Typically the source is expanded in terms of spherical 
harmonics, and then a similar expansion for •0 ret is used 

^ ret = £$£( r '*)iWM) (4) 

where ^m( r ? is found numerically. The individual £m components of , ret in this expansion 
are finite at the location of the particle even though their sum is singular. Then T]^ is 
finite and results from summing gV ' a (tpl^Yi m ) over m. The £-mode expansion of T^ijp) 
was initially determined by a local analysis of the Green's function for an orbit at r Q in 
Schwarzschild coordinates in Ref. [10], 

Inn = (t + 04, + B a + + 0(r 2 ), (5) 

in which it was found that the 0(£~ 2 ) terms yield precisely zero when summed over £. 
Moreover, for circular geodesies in the equatorial plane of the Schwarzschild geometry, the 
regularization parameter C a = and A a and B a also vanish except for their r components. 
The values of A r and B r , first determined by Barack and Ori [10-12], are given below in 
Eqs. (13) and (14). A further term, which we shall denote as D' a , was also introduced in 
Ref. [10-12] and shown there to be zero. It refers to the sum of the 0(£~ 2 ) terms in Eq. (5). 
We comment further about the contribution of D' a towards the end of this section. The 
self-force is ultimately calculated as 




Burko [15] notes in his numerical analysis that the terms in this sum scale as I /I 2 for large 
I, the sum converges as l/£, and it is evident from his results that he computes to at least 
£ = 80 and finds improved convergence with Richardson extrapolation. 



4 



From our perspective, the self-force at a point p' on T is formally given by 



(7) 



evaluated at the source point p'. Formally, the function ip s is defined only in a neighborhood 
of p'; however for calculational purposes, the function may be extended in any smooth man- 
ner throughout the spacetime. While the spherical harmonic components of this extended 
function in the Schwarzschild geometry are not uniquely determined, they still provide a 
convergent expression for ip s for events near p'. Thus, in the Schwarzschild geometry the 
spherical harmonic expansions of ip s and ip ret yield 



and the self-force can be determined by evaluating the vector field 

im 

= V a J2mm-^m)Ytm 



Im 



at the source point p' . Further, with the definitions 



T 



S/ret 



la 



S/ret yr 



(8) 



(9) 



(10) 



the self-force is 



(ii) 



evaluated at r Q . In the above expressions the difference in multipole moments must be taken 
before the summation over t. 

In our approach the regularization parameters are derived from the multipole components 
of V a ^ s evaluated at the source point and are used to control both singular behavior and 
differentiability. In Section V we consider circular orbits of the Schwarzschild geometry at 
radius r n and show that 



lim T% 



£ + 



2y/2D r 



(2£- l)(2£ + 3) 
E l r V 3/2 



(2£-3)(2£- l)(2£ + 3)(2£ + 5) 
where the regularization parameters are independent of £ and given by 

[r (r -3M)]V2 



+ 0(f 6 ). 



A r = -sgn(A)- 



r 2 (V n 



2M) 



B r = 



r Q -3M 
r 4 (r - 2M) 



1/2 



'1/2 



3M)K 



3/2 



2(r - 2M) 



(12) 



(13) 
(14) 
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and 



M(r Q - 2M)F_ 1/2 (r Q - M) (r Q - 4M)F 1/2 
2r 4 (r - 3M) 8r 4 (r - 2M) 

(r -3M)(5r 2 -7r M- 
+ 16r 4 (r - 2M) 2 

-E 1 * has not yet been determined analytically, but the constant V3/2 in Eq. (12) is independent 
of I and is given in Eq. (D23); the £ dependence of the E\ term, and of higher order 
parameters (E*, k > 1), is discussed in Section V. In these expressions F q refers to the 
hypergeometric function 2 Fi [q, \] 1; M/(r Q — 2M)\. The A r and B r terms agree with the 
results of references [10-13] restricted to circular orbits. When summed over all £, the D r 
and terms individually give no contribution to the self-force. This is consistent with the 
results in [10-13], but note the different definition of D r there, which we have referred to 
above as D' a . Our results thus yield the identical self-force to that of Barack and Ori. 

As we shall show in Section IV, in general ip s can be known only approximately. If we 
ignore the D r and terms in the approximation for ifj s , then the approximation for ip R is 
only C 1 . Hence, the D r and E^ terms must be included for ^ K to be a homogeneous solution 
of Eq. (30) as discussed above. Although we have just indicated above that the D T and E* 
terms give no overall contribution to the self-force, we find that understanding the nature 
of these additional terms can be used to speed up dramatically the convergence of the sum 
in Eq. (11). We have used this understanding in obtaining the results of Section E. 

III. THZ NORMAL COORDINATES 

The scalar wave equation takes a simple form when written in a particular coordinate 
system in which the background geometry looks as flat as possible. Consider a geodesic T 
through a background vacuum spacetime geometry g ab . Let 1Z be a representative length 
scale of the background geometry — the smallest of the radius of curvature, the scale of 
inhomogeneities, and the time scale for changes in curvature along Y. A normal coordinate 
system can always be found [17] where, on T, the metric and its first derivatives match the 
Minkowski metric, and the coordinate t measures the proper time. Normal coordinates for a 
geodesic are not unique, and we use particular coordinates which were introduced by Thorne 
and Hartle [7] and extended by Zhang [8] to describe the external multipole moments of a 
vacuum solution of the Einstein equations. In Appendix A, we give a constructive algorithm 
for finding these THZ coordinates for any particular geodesic in a vacuum spacetime. In 
THZ coordinates 

9ab = Vab + H ab 

= Vab + iEab + 3 H ab + 0(p 4 /n 4 ), p/K - 0, (16) 



2r \r - 2M) 
r G -3M 



UM z )F 3/2 3(r -3MY(r + M)F 5/2 



16r 4 (r - 2Mf 



■ (15) 
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with 



2 H ab dx a dx b = 



£ ij x i ^(dt 2 + 5 ki dx k dx l ) 
4 

20 

SijX X^Xfc — p Sij^X 



I ^ a fgpq 13 %x x dt dx 

2 
5 



+ 



21 

5 

21 



dt dx k 



^i^jpq^ fcXP X O Cpqi^Sj X 

O 



dx 1 dx 3 . 



(17) 



and 



3 H ab dx a dx b 



-^S ijk x l x j x k (dt 2 + 5 M dx k dx l ) 

+ ^p^i^x^dt dx k + 0(p 4 /n% dx 1 dx\ 



(18) 



where i] ab is the flat Minkowski metric in the THZ coordinates (t,x,y,z), is the flat 
space Levi-Civita tensor, p 2 = x 2 + y 2 + z 2 and the indices i, j, k, I, p and q are spatial 
and raised and lowered with the three dimensional flat space metric 5^. Note that a term of 
0(p A /TZ A ) is only known to be C 3 in the limit. We call coordinates where H ab matches only 
Eq. (17) second order THZ coordinates; these coordinates are well defined up to the addition 
of arbitrary functions of 0(p 4 /1Z 3 ). Third order THZ coordinates match Eq. (16) through 
the terms in Eq. (18); these are well defined up to the addition of arbitrary functions of 
0(p b /TZ 4 ). The external multipole moments are spatial, symmetric, tracefree tensors and 
are related to the Riemann tensor evaluated on T by 



R 



titj i 



and 



Y2 — „P1f) /O 
£ijk = [^kRtitj] S 



3 r „„_ „ , STF 



pqjt\ 



(19) 

(20) 
(21) 

(22) 



where STF means to take the symmetric, tracefree part with respect to the spatial indices 
i, j and k. Sij and Bij are 0(1/TZ 2 ), and Sijk and Bijk are Oil/lZ 3 ). The dot denotes 
differentiation of the multipole moment with respect to t along T. That all of the above 
external multipole moments are tracefree follows from the assumption that the background 
geometry is a vacuum solution of the Einstein equations. 

The THZ coordinates are a specialization of harmonic coordinates, and it is useful to 
define the "Gothic" form of the metric 



g ab = ^ g ab 



as well as 

H ab = r] ab — g ab 

A coordinate system is harmonic if and only if 



d a H ab = 0. 



(23) 
(24) 

(25) 
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Zhang [8] gives an expansion of g ab for an arbitrary solution of the vacuum Einstein equations 
in THZ coordinates, his equation (3.26). The lower order terms of H ab in this expansion are 



where 



and 



gab = ^ab + ^ab + (p 4 /^ 4 ) 



Tab 



ab 



(26) 



yH tk 



~ 2 Es<i j x X J 



2 10 

6? JSn^ X TiX ' I 

3 q p 21 



j *^ ^ p> *^ P 



2- 
—i 

5 



5 

21 



x 



(t e^B qk x p x k - \e p ^&\x pP 2 
o 



(27) 



H tk 
gij 



^6 ^BqijXpX '' 3? 

o( P A /n% 



(28) 



At linear order in H ab , the metric perturbation H ab is the trace reversed version of H ab , 



ab 



H a b — H ab — -g ab H c c , 



(29) 



and Eqs. (16)-(18) are precisely the terms up to 0(p 4 /7£ 4 ) which correspond to Zhang's [8] 
expansion. 



IV. GREEN'S FUNCTIONS FOR A SCALAR FIELD 

The scalar field equation 

V 2 -0 = -Ang (30) 
is formally solved in terms of a Green's function, 

v 2 GW) = -(-g)- 1/2 ^(x; - x ;,), (31) 

where p' represents a source point on T, and p a nearby field point. The source function for 
a point charge moving along a worldline T, described by p'(r), is 

g(p) = qj{-g)- 1/2 5\p-p\r))dr, (32) 

where r is the proper time along the worldline of the particle with scalar charge q. The 
scalar field of this particle is 

^(p) = Airq [ G[p,p'(r)]dr. (33) 
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DeWitt and Brehme [3] analyze scalar-field self-force effects by using the Hadamard 
expansion of the Green's function. An important quantity is Synge's [9] "world function" 
a(p, p') which is half of the square of the distance along a geodesic between two nearby points 
p and p' . The usual symmetric scalar field Green's function is derived from the Hadamard 
form to be [3] 

G sym (p, P ') = ±- [u(p,p')5(a) - v (p,p')Q(-a)] (34) 

where u(p,p') and v(p,p') are bi-scalars described by DeWitt and Brehme. The 0(— a) 
guarantees that only when p and p' are timelike related is there a contribution from v(p,p'). 
The retarded and advanced Green's functions are 

G^(p,p') = 2&[i:(p),p'}G s ^(p,p') 

G^(p,p') = 2Q\p> ^{p)\G s ^{p,p') (35) 

where 0[E(p),p'] = 1 — 0[p', £(p)] equals 1 if p' is in the past of a spacelike hypersurface 
S(p) that intersects p, and is zero otherwise. The terms in a Green's function containing u 
and v are commonly referred to as the "direct" and "tail" parts, respectively. 
A second symmetric Green's function [1] 

G\p,p') = ^ [u(p,p')5(a) +v(p,p')&(a)} (36) 

precisely identifies the part of ip which is not responsible for the self-force. In particular, the 
regular remainder 

is a homogeneous solution of the the field equation (30) and completely provides the self-force 
when put on the right hand side of Eq. (1) [1]. 



A. Approximation for ip' 



In this section approximate expansions are derived for G s and ip s . For a vacuum spacetime 
(R a b = 0) which is nearly flat Thorne and Kovacs [18] show that 

u(p,p') = l + 0(p'/n'), (38) 

their equations (39) and (40), and they evaluate the direct part of the retarded Green's 
function to be 

±-u(p,p') 5 TCt {a(p,p>)} = ( 1 + °^l /n4) ) S(t p -t rct ), (39) 

V / Trot 

where the dot denotes a derivative with respect to t p > 

We now express & TCt in terms of the THZ coordinates to obtain Eq. (47) below. When 
the source point p' is on T, a(p,p') is particularly easy to evaluate in THZ coordinates for p 
close to p'. Synge's [9] "world function" a(p,p') is shown by Thorne and Kovacs [18] to be 

a(p,p') = l -x a x h (^ ab + J H ab dX^j + 0(p 6 /n 4 ), (40) 
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their equations (37) and (38), where the THZ coordinates of p' are (t p /, 0,0,0), x a is the 
coordinate of the field point p, and the coordinates of the path of integration C are given 
by ( a (\) — M 2 -" ~~ V^") w ith A running from to 1. We closely follow the analysis in [18], 
while using THZ coordinates, and only work through lower orders in p/TZ. 

Given H ab = 2H ab + zH ab from Eqs. (17) and (18), the integral of a component of H ab 
along C is straightforward. For example, 

J H tt d\ = - ^(svCC + l^kCVC^ d\ + o( P A /n A ) 

= -U, r r'.r' - i^xV!* + 0(p 4 /n 4 ), (41) 
The other components give similar results. If we define 

H ab = ^H ab dX, (42) 

then Synge's world function is 



<j(p,p') = -x a x b rj ab + -(t p - t p ,) 2 H tt + (t p - t p/ )x l H tt + -x^Hij + 0(p 6 /TZ 4 ), 

= ~(i - H tt ) [(v -t P + x l n lt ) 2 - ./•'./•'!//„ + n t ,)/(i - n tt ) 

+ 0(p 6 /n 4 ). (43) 

The second equality depends upon the facts that H ab = 0(p 2 /1Z 2 ) and that \t p ' — t p \ = O(p) 
near the null cone. 

With the source point on T, 

xVfao- + Hij) = P 2 (l + Hu) + 0(p 6 /TZ 4 ) (44) 

where 

p 2 = x l x j rjij. (45) 

The result in Eq. (44) depends upon the detailed nature of 2 Hy and sH^ in Eqs. (17) and 
(18) as well as upon the definition of Hij in Eq. (42). 

After the substitution of Eq. (44) into Eq. (43), factorization of a yields 

<r(p,p') = —(l-HttWp'-tp + x'Hit-pil + Htt)} 

x [V -t p + xm it + p(l + Hu)] + 0(p 6 /^ 4 ). (46) 

At the retarded time, p' is on the past null cone emanating from p, where cr(p,p') = 0, 
and it follows that the first of the factors in square brackets is ~ p and the second must 
be p' 1 x 0(p 6 /1Z 4 ) = 0(p 5 /1Z 4 ) to cancel the order term in Eq. (46) and have a(p,p') 
vanish precisely. Thus, differentiation of Eq. (46) with respect to t p > and evaluation at the 
retarded time yields an expression which is dominated by the part which results from the 
differentiation of the second term in square brackets, 

~da(p,p') 



dtp/ 



= - w«)[V -t P + x l H tt - p(l + n tt )] + 0(p 6 /n 5 ) 

ret 

= -l(l-H tt )l-2p(l + H tt ) + 0(p 5 /n 4 )} 

= p[l + 0(p 4 /n 4 )}; (47) 
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the first equality follows from taking the derivative of the second term in square brackets in 
Eq. (46) with respect to t p r, the second equality from evaluating at the retarded time, and 
the third equality follows from Eqs. (41) and (42). 

The direct part of the retarded Green's function in Eq. (39) is now 

-^-u(p,p') S TCt [a(p,p')] =(4npy 1 5[t p/ -t p + x l Uit 

+ p(l + Hu) + 0(p 5 /n 4 )] [1 + 0{p'/n% (48) 

DeWitt and Brehme show that in general 

v(p,p') = —^R(p') + 0(p/n 3 ), p^T, (49) 
but in vacuum, where R = from the Einstein equations, 

v(p,p') = o( P 2 /n 4 ). (so) 

This follows from Eq. (2.9) with substitutions from Eqs. (2.14), (2.15), (1.76) and (1.10) 
of Ref. [3]. The dominant contribution to ip s from the v(p,p') term is 0{p i /TZ A ) in the 
coincidence limit, p — > V . 

All together then, with Eqs. (48) and (49), substituted into Eq. (36) and an integration 
over the worldline, 

^ s = q/p + 0(p 3 /n 4 ). (51) 

We note that the third order THZ coordinates are only well defined up to the addition of 
a term of 0(p 5 /TZ A ). Such an addition would change 1/p by the sum of a term that is 
0(p 3 /7?. 4 ), and would be consistent with the order term of Eq. (51). The differentiability 
of the order term is of interest, and a term of 0(p 3 /7?. 4 ) is C 2 in the limit that p — > 0. 
In light of the fact that ip R is a homogeneous solution, Eq. (51) clarifies the relationship 
between the accuracy of an approximation for ip s and the differentiability of the subsequent 
approximation for ip R = ip ret — ip s , and the self-force d a ip R . Specifically, if the approximation 
for ?p s is in error by a C n function, then the approximation for ?/> R is no more differentiate 
than C n and the approximation for d a i[) K is no more differentiable than C n ~ l . 



B. Intuitive understanding for ip° 

Before continuing, it is instructive to provide an elementary, direct explanation of Eq. (51) 
by taking full advantage of the features of the THZ coordinates. The scalar wave operator 
in THZ coordinates, is 

v^V a V ^ = d a (77^) - d a (H a %ij) (52) 

or 

v^V a V ^ = v ab dad b *p - H ij did^ - 2H u d {i d t) iP - H u d t d t ^. (53) 

Direct substitution into Eq. (53) shows how well q/p approximates ?/> s . If ip is replaced 
by q/p on the right hand side, then the first term gives a 5-function, the third and fourth 
terms vanish because p is independent of t, and in the second term 2 H % i has no contribution 



11 



because of the details given in Eq. (27), and the remainder of H l i yields a term that scales 
as 0(p/K 4 ). Thus, 

^V a V a {q/p) = -Anq5\x l ) + 0{p/TZ A ), p/TZ - 0. (54) 

From consideration of solutions of Laplace's equation in flat spacetime, it follows that a C 2 
correction to q/p, of 0(p 3 /TZ 4 ), would remove the remainder on the right hand side. We 
conclude that ip s = q/p + 0(p 3 /TZ 4 ) is an inhomogeneous solution of the scalar field wave 
equation. And the error in the approximation of ip s by q/p is C 2 . 



V. REGULARIZATION PARAMETERS FOR A CIRCULAR ORBIT OF THE 
SCHWARZSCHILD GEOMETRY 

In Appendix B we give the detailed functional relationship between the THZ coordinates 
(t, x, y, z) and the Schwarzschild coordinates (t s , r, 9, 0) for an orbit described by r = r Q with 
6 = 7r/2 and = Vlt s . 

As seen in Section IV, an approximation to ip s is 

^ = q/p + 0(p 3 /n 4 ). (55) 

The regularization parameters result from evaluating the multipole components of q/p at 
the location of the source. The use in this manner of q/p, in lieu of ip s itself, is justified 
because the error in the approximation to ip s , being 0(p 3 /7?. 4 ), gives no contribution to 
V a ip s as x — > 0. 

To aid in the multipole expansion we rotate the usual Schwarzschild coordinates to move 
the coordinate location of the particle from the equatorial plane to a location where sin 6 = 
for a specific t s , following the approach of Barack and Ori as described in [11]. Thus, we 
define new angles 9 and $ in terms of the usual Schwarzschild angles by 

sin 9 cos(0 — fit s ) = cos 6 

sin 9 sin (0 — Qt s ) = sin© cos $ 

cos 9 = sin© sin $. (56) 

A coordinate rotation maps each Yg m (9, (j>) into a linear combination of the Y^ m /(0, $) which 
preserves the index £, while m! runs over —£...£. Thus, the £ component of the self-force, 
after summation over m, is invariant under the coordinate rotation. 

To obtain the regularization parameters: first we expand d r (q/p) into a sum of spherical 
harmonic components whose amplitudes depend upon r. Then we take the limit r — > r a . 
Finally only the m = components contribute to the self-force at = because Yf m (0, $) = 
for m 7^ 0. Thus, the regularization parameters of Eq. (12) are the (£,m = 0) spherical 
harmonic components of d r (q/p) evaluated at r . 

In this section e is a formal parameter which is to be set to unity at the end of a calculation; 
a term containing a factor of e n is 0(p n ). We use it to help identify the behavior of certain 
terms in the coincidence limit, x l — > 0. We have used Maple and GRTensor extensively 
to obtain the results reported below. 

A lengthy expression for p 2 , for a circular orbit in the Schwarzschild geometry, may be 
derived from the analysis of Appendix B. The 0(e 2 ) part of p 2 is 
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where 



and 



x = i- 



(r - r Q ), 
M sin 2 $ 



(58) 



(59) 



r G -2M 

The dependence of p 2 on the Schwarzschild coordinates may be written solely in terms of 
A, x and P by use of Eq. (57) to remove cos0, then Eqs. (58) and (59) remove r and sin$ 
respectively. We formally expand d r (l/p) in powers of e to obtain 



d r (l/p) = 



p 3 (r G -2M) pr 
+ e °0(A/p,A 2 /p 2 ) + 



3M 



+ 0(A 2 /p 2 ,A 4 /p 4 ) 



2 X (r Q - 2M) 

M(r -2M) (r -M)(r -4M) 



2(r G - 3M) 



+ 



(r G - 3M)(5r G 2 - 7r M - 14M 2 ) 

16x 2 (r G - 2M) 2 
3(r Q - 3M) 2 (r G + M) 



+ 0(A 2 /r o ,A 4 /p 2 rc) 



16x 3 (r G - 2M) 2 



8 X (r - 2M) 



+ 0(p 2 ). 



(60) 



We consider the multipole expansion, in the limit that A — > 0, of the m = part of 
each coefficient of e n for n = —2 to 1. A convenient method to find the m = component 
in one of the following terms involves integrating over the angle $ using details described 
in Appendix C. The expansion of the 6 dependence in terms of Legendre polynomials is 
described in Appendix D. 

First term The A — > limit of the expansion of the e~ 2 term in Eq. (60) is 
A% p3(r - 2M) 



1/2 



lim — . 

A^o \r - 2M 
r -2M) ,■ 



r n A 5 



r -2M 



1/2 r 



,A 2 



r D -2M 



+ 



2r 2 (r -2M) 
r G -3M 



X(l - cos 6) 



S gn(A) (r „-3M) f./ <+ l) Pt(cM9)t 
2 (r -2M-M S in 2 *)f-;V 2/ " 



-3/2 



(61) 



where e has been set equal to 1 on the right hand side here and below, and where the second 
equality follows from Eq. (D12) with the substitution 



5 2 = A 2 (r G - 3M)/[2r (r - 2M) 2 X ]. 



(62) 



Integrating over <l> and dividing by 2tt (denoted by the angle brackets () here and in Ap- 
pendix C) via Eq. (C7), to find the m = contribution, results in 



I lin o\ ~p\r -2M)j 



e- 2 r D A \ 



= -sgn(A) 



[r (r -3M)] 1 / 2 
r 2 (r - 2M) 



E(^ + 0^(cos0). 



(63) 
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In the coincidence limit P^(cosO) = 1 and a term in this sum is then 

-^+^sgn(A) Mr °- 3M)1 ' /2 
{ e+ 2 ) Sen[a > r,'(r„-2M) 

which determines the A r term in Eq. (12) as given in Eq. (13). 
Second term The A — > limit of the e" 1 term in Eq. (60) is 



lim — — 

A^o pr Q 



1 - 



r G -3M 
2 X (r Q - 2M) 



r G -3M 



2r 4 (r -2M- M sin 2 $)(1 - cos 9) 



1/2 



1 - 



r Q -3M 
2 X (r G - 2M) 



r G -3M 
2r 4 (r - 2M) 



1/2 



1 - cos 9) 



-1/2 



r D -3M 



x i/2 2x 3/2 (r Q -2M) 



(64) 



(65) 



Integrating over <£> results in hypergeometric functions as shown in Eq. (C3), and the ex- 
pansion of (1 — cos9)~ 1 / 2 in terms of the P£(cos9) is given in Eq. (D7) and results in 



lim ( — — 



r Q -3M 
X {r -2M)\ 

1/2 



r Q -3M 
2r 4 (r - 2M) 



'1/2 



(r Q - 3M)F 3/2 
2(r Q - 2M) 



(66) 



In the coincidence limit P^(cos9) = 1 and a term in this sum is then 



B r = 



r G - 3M 
[r Q \r Q -2M)\ 



1/2 



'1/2 



(r Q - 3M)F; 



3/2 



r G - 2M 



which is the B r term in Eq. (12) as given in Eq. (14). 



(67) 



Third term The 0(e°) term in Eq. (60) is zero in the limit that A — > for nonzero 9, 
and gives no contribution to the sum in Eq. (12) as follows from Eq. (D7). 



Last term For the last, e 1 , term in Eq. (60) we consider 



lim p 



2r 2 (r - 2M) 



r G -3M 



1/2 



X 1/2 (l-cos9) 1 



/2 



(68) 



After the expansion of (1 — cos 9) 1 / 2 described in Eq. (D16) and the integration over $ with 
Eq. (C3), the multipole expansion of the p terms in Eq. (60) gives 



2r 2 (r -2M) 
r -3M 



1/2 



M(r G - 2M)F_ 1/2 (r D - M) (r Q - 4M)F, 



1/2 



2r 4 (r - 3M) 



8r 4 (r - 2M) 



(r - 3M)(5r 2 - 7r M - 14M 2 )F 3/2 3(r G - 3M) 2 (r + M)F 5/2 



x 



E 



16r 4 (r - 2M) 2 
-2V/2PKCOS9) 



16r 4 (r - 2M) 2 



(2£- l)(2£ + 3)' 



(69) 
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In the coincidence limit P^(cos6) = 1, and a term in this sum is then 



-2\f2 

Dr (2£-1)(2£ + 3Y (70) 
which defines D r , as in Eq. (15), and is a term in Eq. (12). 

The remainder These derivations of the regularization parameters reveal a pattern 
for the ^-dependence of higher order parameters, even if the overall scale of the parameter 
remains unknown. 

The successive terms in Eq. (60) provide increasingly accurate approximations of d r ip s 
and also increasingly accurate approximations of d r ifj R = d r ifj ret — d r ip s . In principle the 
terms through O(e ) are sufficient to calculate the radial component of the self-force; this 
is effectively the level of approximation described in References [10-13] and implemented in 
Ref. [15]. With the inclusion of 0(e°) terms the approximation for d r ip R is C°, the remainder 
terms scale as £~ 2 for large i and their sum converges. But when the approximation of d r ip is 
improved by the addition of the 0(e) terms, the resulting approximation of d r ip R is then C 1 , 
and we see below that the remainder terms scale as £~ A resulting in a more rapid convergence 
of the sum for the self-force. 

As the approximation of d r ip s is improved by successive terms of greater differentiabil- 
ity, the resulting approximation of d r ip R is not only more differentiate but also leads to 
increasingly rapid convergence of the sum for the self- force. 

We can anticipate the details of how this occurs. From the descriptions of the THZ 
coordinates in Appendix B and of ip s in terms of THZ coordinates in Section IV, we expect 
that the 0(e 2 ) term in Eq. (60) is C 1 . A more accurate approximation to ip s could be 
provided by the modification 

p 2 - p 2 + X N X N , (71) 

where the components of A at = 0(l/1Z n ~ 2 ) are not functions of the coordinates but depend 
only upon the orbit and the coordinate location of the particle. Here we borrow the notation 
of the analysis of STF tensors [19] where AMs a multi-index that represents n spatial indices 
. . . i n , however while Atv is symmetric it is not necessarily tracefree. Also X N represents 
X ll X 12 . . . X ln where the X 1 represents one of X, y or Z defined in Appendix B. Now, p 2 
is already determined through 0(e 5 ), so that n is necessarily greater than or equal to 6 for 
an improved approximation. And while we do not know the actual value of A at, we assume 
that such a Aa? exists that provides an improved approximation to ip s . 
Such a correction to p 2 ultimately results in the addition of 



8 r 



X]\[X^ 



(72) 



to d r (l/p) in Eq. (60), which is of 0(e 2 ) and C 1 if n = 6. This result is consistent with 
Eq. (51) and with the discussion following Eq. (54) in Section IV above. With this improved 
approximation the resulting error term in Eq. (60) would be 0(e 3 ) and C 2 . 

Finding higher order THZ coordinates might not be the best way to correct the approx- 
imation for -^ s , but any such correction involving an expansion about T would necessarily 
take the generic form of a homogeneous polynomial in X, y and Z (or equivalently in x l ) 
divided by p raised to some integral power. Thus to find higher order corrections to ^ s , we 
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are led to consider the multipole expansion of a term such as in Eq. (72), for n > 6 and to 
determine the nature of the regularization parameters that would result. 

First X, y and Z are replaced by their definitions in terms of the usual Schwarzschild 
coordinates in Eqs. (B1)-(B3). Then, the angles are changed to 6 and $ via Eq. (56). And 
finally all coordinate dependence is written in terms of A, x an d p as described above in 
Eq. (60). The result is a sum of terms each of which is of 0(e™ -4 ), for n > 6 and whose 
coordinate dependence is contained in the functions A (which depends only upon r), x 
(which depends only upon $) and p. As x — > both A and p are 0(e) while x — 0(1), so 
that each term must include a factor A 9 /? 5 for integers q and p where q + p = n — 4 > 2. In 
fact, careful analysis of the above substitutions shows that q > 0. All of the 6 dependence 
resides in pP. 

The Legendre polynomial expansion of p p , for odd p, is discussed at length in Appendix 
D. There we show that when p — — 1 (respectively, p < —1), a consequence of Eq. (D7) 
[respectively, Eq. (Dll)] is that the expansion coefficients scale as a constant (respectively, 
diverge as (r — r ) p ) when r — > r . The coefficients of the expansion A^/F then scale as 
(r — r Q ) q (respectively, (r — r Q ) q+p ). In both cases the power is an integer > 2. In the 
limit that r — > r Q these terms approach zero and give no contribution to the regularization 
parameters. The case of p even and negative gives a similar result, but is not discussed in the 
Appendix. We conclude that if p < then no contribution to the regularization parameters 
results. 

If p > and q > then the Legendre polynomial expansion of p p is well-behaved and 
finite as r — * r G , but the product A 9 /^ vanishes in the limit r — > r G and gives no contribution 
to the regularization parameters. 

If q = and p > 2 and is even, then the G dependence is in the form of a polynomial in 
cos0 which has an expansion in terms of the Legendre polynomials only up to P P (Q), and 
because p = when r = r Q and 6 = the sum of this finite number of terms is zero and 
gives no contribution to the regularization parameters. This case always results when the 
improvement to d r ip s is 0(e n ) for n being even. 

The only remaining case is q = and p > 2 being a positive odd integer. In the limit 
that r — > 0, ff oc (1 — cos©) p / 2 . The Legendre polynomial expansion of this function is 
discussed in detail in Appendix D. We see in Eq. (D22) that, for k a positive integer 

oo 

(1 - cos0) fc+1 / 2 = ^^ +1/2 P,(cos0) (73) 
e=o 

where 

A k e +1/2 = (2£ + l)V k+1/2 / [(2£ -2k-l)(2£-2k + l)... (2t + 2k + \)(2t + 2k + 3)] , (74) 

for a constant Vk+1/2 given in Eq. (D23). When 9 = 0, P£(cos6) = 1 and such terms 
do contribute additional regularization parameters in the mode sum representation of the 
self- force. This case always results when the improvement to d r vp s is 0(e n ) for n being odd. 

We now see that every other higher order correction to ip s provides an additional regu- 
larization parameter. The i dependence is necessarily of the form 

E k a 4 +1/2 (75) 

where E% is independent of £ but still undetermined. The first term of this sort, for k = 1 
is included in Eq. (12). It important to note that for each value of k the sum of these terms 
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TABLE I: The fitted parameters of F\f 
k £* part of ^ r R 

1 1. 80504(4) x 10" 4 2.78201(7) x 10~ 9 

2 -1.000(3) x 10~ 4 6.90(2) x 10~ 12 

3 4. 3(3) x 10~ 5 3.1(2) x 10~ 14 

4 -5. 6(5) x 10~ 5 7.7(6) x lO" 16 

from £ = to infinity is necessarily zero, and need not be included in the self-force analysis. 
However we see in the next section that including these additional coefficients dramatically 
speeds up the convergence of the self-force sum Eq. (11). 



VI. APPLICATION 

In this section we apply the formalism developed above to determine the self-force JF R 
for a scalar charge in orbit at r = 10M about a Schwarzschild black hole. In the numerical 
work we use units where M — 1 and q — 1. Appendix E describes the practical details for 
numerically integrating the scalar wave equation to determine the J~lf"- To compute the 
self- force, the A r , B r and D r terms must first be removed from as in Eqs. (11) and 
(12) — a process which determines residuals which fall off as £~ 4 for large £. Removing the 
contribution of each successive improves the falloff of the residuals by an additional two 
powers of £. 

If we have Tlf 1 for all £, summing the residuals after removing A r , B r and D r would give 
us the self-force. With evaluated only for finite £, we must make some attempt to obtain 
the higher £ contributions. To do this we will numerically determine the coefficients by 
fitting the residuals, considered as a function of £, with a linear combination of terms whose 
£ dependence is given by the E\ term in Eq. (12) for k — 1 and by Eq. (75) for integers 
k>l. 

A comparison of the integration results, J-Jf, for different values of an accuracy parameter 
in the numerical routine, revealed that a systematic effect remained in our best data for T\f . 
In order to avoid fitting to that systematic effect, we chose to add to our data a small random 
component which was capable of swamping any trace of the systematic effect, and which 
allowed us to have precise control of the error in the T\f . This random component also 
provided the opportunity to use Monte Carlo analysis to determine the statistical significance 
of our result of the self-force. 

In fitting for the we avoided small values of £, which may contain significant physical 
information not associated with the large £ falloff of J-Jf 1 - Thus we fitted the residues for 
£ from 13 to 40 while determining from 1 to 5 of the coefficients. Fig. 1 summarizes 
the results of this numerical analysis. The curve labeled Fn is as a function of £. The 
curves A, B and D show — Tf r where J-f r successively includes the contribution from 
the regularization parameters A r , B r and D r . The E 1 to E 5 curves show the residuals 
after numerically fitting from 1 to 5 of the E^ coefficients and removing their contributions 
successively. 



We actually used a singular value decomposition from Ref. [20] to fit the residuals. It 
provided an independent estimate of the uncertainty of the E* 's, which is entirely compatible 
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FIG. 1: The upper portion of the figure displays as a function of £, along with the result of 
it being regularized by A r , B r , and D r . The lower portion displays the residual after a numerical 
fit of from 1 to 5 additional parameters, E\ . . . E^. A point where the data on a particular curve 
changes sign from being negative to positive is labeled with +, from positive to negative by o. 

with the Monte Carlo results. This represented a valuable, overall consistency check of our 
analysis. The E^ coefficients which result from a fit of four coefficients are given in Table 
I along with their uncertainties in brackets. After fitting four coefficients and removing 
their contribution, we obtained an RMS residual of 2.8 x 10~ 14 over the fitting range, which 
is completely determined by the size of the random component we had introduced to the 
original J^f to swamp the systematic effect. Four coefficients evidently fit the data down 
to the noise. It was clear that fitting a fifth coefficient, or more, did not improve the quality 
of the fit. 

The self-force JF^ = 1.37844828(2) x 10~ 5 was obtained by summing, over the range of our 
data, T^ 1 with the A r , B r and D r terms removed as in Eqs. (11) and (12). The remainder 
of the sum to £ = oo was approximated by the contributions of the ££, E%, ... sums from 
41 to oo, once the coefficients had been determined. The uncertainty was obtained from 
the Monte Carlo simulation. The table also shows the individual contribution of each 
to the self-force Tf- as well as the amount of the uncertainty in Ff- which is attributable to 
that E*. Without including the effects of the £* tails, we would have found the sum out 
to £ = 40 to be 1.37817 x 1CT 5 . Fitting the higher order terms has allowed us to increase 
dramatically the effective convergence to our final result. 

Our result is consistent with Fig. (4A) of Burko's analysis[15]. With the A r and B r terms 
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removed from J-^f, he effectively calculates the total self- force by summing data points on 
the equivalent of our curve B out to a large enough value of £ that convergence is obtained 
while using Richardson extrapolation. 

A future manuscript will apply our methods to the investigation of physical questions. 

VII. DISCUSSION 

In a previous paper [1], we outlined our method for computing the self-force. This hinged 
on realizing that ip R = -?/> ret — ?/> s is a homogeneous solution of the field equation and we proved 
that it gave the same result as methods based on using the tail part of the retarded Green's 
function. As a consequence of this we have obtained the same regularization parameters as 
all previous authors [10-13] using a regularization procedure based on mode sum expansions. 
Exact computation of ip R would yield a homogeneous solution of the field equation. Under 
interesting physical circumstances, we anticipate that ip R should consequently have a high 
level of differentiability [31]. The level of differentiability of an approximation for ^ R is 
limited by the accuracy of the approximation for vp s . To improve the level of differentiability 
of our approximation for ip R beyond that of ?/> tai1 , we have thus been led to explore higher 
order approximations to ip s . 

Following earlier work [21], we have used THZ coordinates to obtain the simple approx- 
imation ip s q/p. Our key analytical result is the expansion in Eq. (60) which is based 
on this approximation. The regularization parameters are derived from the mode sum rep- 
resentation of each term in Eq. (60). The parameters A r and B r come from the first two 
terms. The parameter C r is seen to be zero, directly from the third term. All these results 
are consistent with previous work by others [10-13]. Our D r parameter comes from the 
fourth term in Eq. (60), which we compute analytically and for which we find the specific £ 
dependence of l/(2£ — l)(2£ + 3). Direct inspection of this term shows that the sum goes to 
zero in the coincidence limit, so it does not contribute to the self-force. Nevertheless, it is 
recognition of the large ^-behavior of the mode sum expansion for this, and similar higher 
order terms characterized by the E*, which leads to dramatically improved convergence in 
the mode summation. Understanding of the specific nature of the £ dependence in the mode 
sum representation of the higher order terms was obtained by an analysis of general methods 
for improving the approximation to ip s . Our numerical application of this scheme amply 
illustrates the benefits of estimating the parameters in order to accelerate convergence. 

In principle, neither the use of ip R instead of -?/> tai1 , nor the specific use of THZ coordinates 
are intrinsically necessary in the computation of the self-force. Indeed, many authors have 
used neither, and have yet obtained analytical results for the regularization parameters, 
and/or numerical results for the self-force [10-14, 22-24]. It is clear that other methods 
might also be used to calculate D r as well as the general £ dependence of our terms. What 
we believe is important is that the form of these higher order terms has been determined, and 
that the inclusion of these terms has a dramatic impact on the effectiveness and accuracy 
of numerical work. 

Unquestionably, the intensity of work in this area is paying off, and as efficient computa- 
tional techniques are recognized and implemented, a greater volume of results will become 
available. This will be especially true in relation to the long awaited detection of gravita- 
tional waves from binary inspiral sources, represented by a small compact object in orbit 
about a comparatively large black hole. 
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APPENDIX A: THE DETERMINATION OF THZ COORDINATES 

A particular THZ coordinate system (t, x, y, z) is associated with any given geodesic 
r(r) of a vacuum spacetime and is a harmonic coordinate system as well as being "locally 
inertial and Cartesian." By Zhang's [8] definition a "locally inertial and Cartesian" (LIC) 
coordinate system has the spatial origin x l = on the worldline T, and has the metric being 
expandable about T in powers of p in a particular form which we describe as g a b = f] a b + fP x 
homogeneous polynomials in x l of degree q, for non-negative integers p and q with p + q > 2. 

The defining features of nth order THZ coordinates are that 

(i) On T: t measures the proper time along the geodesic, the spatial coordinates x, y and 
z are all zero, g a b\r = Vab and all of the first derivatives of g a b vanish. 

(ii) At linear, stationary order (c/. Ref. [8]) H % i = 0(x n+1 ). 

(Hi) The coordinates satisfy the harmonic gauge condition d a Q ab = 0(x n ). 

To find THZ coordinates associated with a particular geodesic, it is easiest to satisfy these 
conditions in order. 

Given a general set of coordinates Y A and a particular point p a new set of coordinates 
X a may be defined by 



X a = A a + B a A (Y A -Y p 
1 
2 



+^B a A T A BC {Y B -Y p B ){Y c -Y p c ) + 0[{Y -Y p f] Y A — > Y p A , (Al) 



where T A BC are the usual Christoffel symbols, the A a are arbitrary constants, and the B a A 
are also arbitrary constants restricted by the condition that B a A, considered a matrix, be 
invertible. Weinberg [25] shows that 

y y q Y a dY s 

= r ] ab + 0[(Y-Y p )% Y A — > Y A , (A2) 

a ab 

^- c =0[(Y-Y p )], Y A — > Y A . (A3) 

A specific choice for the A a and B a a result in the coordinates of p being X p = r p 5". And 
condition (i) is satisfied by repeating this construction along T while parallel propagating 
the coordinate basis. 



so that 
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Thus, the coordinates that satisfy condition (i) are denoted X a . And 

g ab = r] ab — H ab 

= rj ab - H ab ij X i X j + 0(X 3 /n 3 ), X i -> 0, (A4) 

here X in the order term can refer to any of the spatial coordinates, and the 

1 d 2 H ab 



(A5) 



J^ a o _ 

ij ~ 2 dXtdXi 

are functions only of t. 

To satisfy conditions (ii) and (Hi) we use a gauge transformation of the form 

X (new) = ^(old) + C° (A6) 

where changes in geometrical objects are calculated only through linear order in ( a . In this 
application, we are interested in the vicinity of T and accordingly let 

C = CijkX i X i X k , (A7) 

where the C a ijk are functions only of t to be determined below. Thus 

xl cw) = X a + ( a ijk X*XiX k , X l ^0 (A8) 

ab 



changes H ab to 



#new = + d% b + d b C - V ab d c C + 0(C 2 ). (A9) 



or 



l -H^ = l -H^ + ( aMj + C bai3 ~ V ab Ck kij - (A10) 

We find the necessary gauge transformation in two steps. To satisfy condition (ii) we 
require that the gauge transformation obey 

dX J + &{< - ^d c ( c = -H% d + 0(X 3 ), (All) 

which implies that 

Cijkl + (jikl - Vij( m mkl = —^H°j^. (A12) 

We use the decomposition of STF tensors of reference [26] and let 

3 5 1 - 

Qjki = A(ij k i) + -eip(jB p k i) + -8i(jCki) — -Hi P p<y j5ki), (A13) 

where (...) implies taking the STF part of the enclosed indices and Aijki, B pk i and Cm are 
STF tensors on all of their indices. The coefficient —1/5 in the last term is chosen to make 
the trace of Eq. (A13) agree with the trace of Eq. (A12). The remainder of the solution for 

Cijkl IS 

Ajki = -g-ff(ijfci> + ^{ijH P kt) P + -^H p p(ij 5 k i) - 7^5(ij5jk)H pq Pq) (A14) 
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Bpki = -^ij( P H\i) J + —5( k ie p ) 1J H iq q j (A15) 

and 

1 1 - 

Cm = -H\ k i + — (H P jk P + H p kjp + 5 jk H pq pg ) , (A16) 

but only if Hy k i satisfies the auxiliary condition that 

Hku — H l k u — H\ ki + 5 k iH pq pq = 0. (A17) 

This condition is automatically satisfied when H ab is derived from a vacuum solution of the 
Einstein equations, cf. Eq. (35.64) of Ref. [17]. 

The spatial component of condition (Hi) is now satisfied through 0(x) because H' tj van- 
ishes. The time component of condition (Hi) will be satisfied through O(x) also if 

d b d b C l = -d b H% d , (A18) 

thus 

d b d b C = -( t ijk X i X j X k + 6( u ik X k = -2H u ik X k + 0(X 2 ) X i -> 0, (A19) 

or 

An elementary solution for C^ijk is 



C\k = -\s u lk . (A20) 



Cijk = --H tp p (i5j k ). (A21) 

The combination of Eqs. (A13) and (A21) provides the necessary gauge transformation 
that results in satisfaction of conditions (ii) and (Hi) for second order THZ coordinates. 

The third order coordinates may be found following a similar procedure where the gauge 
transformation is of the form 

x a {new) ^x a {oiA) +Ci j kiX i X^X k X l , X^O. (A22) 

At the fourth and higher orders, where terms in the metric expansion involving two 
derivatives of the Riemann tensor are of the same order as terms quadratic in the Riemann 
tensor, the presence of nonlinearities complicate the construction of THZ coordinates. We 
have not needed detailed knowledge about these higher order THZ coordinates throughout 
this paper. 



APPENDIX B: THZ COORDINATES FOR A CIRCULAR ORBIT OF 
SCHWARZSCHILD 

To find the THZ coordinates for a circular orbit in the Schwarzschild geometry it is 
convenient to take advantage of the spherical symmetry of the background geometry rather 
than to follow the general procedure described in the preceding section. The orbit T, given 
by = fit where f2 = a/ Mj r Q 3 is the orbital frequency at Schwarzschild radius r Q , is tangent 
to a Killing vector field £ a = d/dt s + £ld/d(j). We choose three helping functions X, y and Z 
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which are Lie derived by the Killing vector, C^X = C{y = C^Z = 0, and their gradients are 
spatial and orthogonal to the 4-velocity of the geodesic and to each other when evaluated 
on T. These helping functions are 



r — r D 



and 



X 

y = r Q sin 9 sin(0 — VLt 

Z = r Q cos(6>). 



(l-2M/r Q ) 1 /2' 

r G - 2M^ 1/2 



r G - 3M 



(Bl) 
(B2) 
(B3) 



x 



Two more useful functions are 
[r sin 9 cos(0 — Ot s ) — r Q ] 



+ 
+ 
+ 



'^-' 2 -v ,9 
h 3^ 

2 y 



(1 - 2M/r Q ) 1 /2 
M 

r G 2 (l -2M/r )V2 
MA 1 

2r 3 (r -2M)(r -3M) 
M 

r 5 (l-2M/r )V2( ro _3M) 



r G - 3M 



\-M 2 X 2 + ^ 2 (r - 3M)(3r G - 8M) + 3Z 2 (r G - 2M) 2 ] 



MX 



x 2 y 2 

+ — f-(28r G 2 - 114r G M + 123M 2 ) + , 



(r Q 2 - r Q M + 3M 2 ) 
8(r Q - 2M) 

I4r 2 - 48r G M + 33M 2 ) 



+ 



My 4 



4(r G - 2M) 



(3r Q 3 - 74r G 2 M + 337r Q M 2 - 430M 3 ) 
(3r G + 22M) 



56(r G - 2M) 2 

M 2 y 2 Z 2 (7r - 18M) MZ 4 



56 



(B4) 



and 



y = r sin#sin(0 — Qt 

Mxy 

+ 



r - 2M 



r Q - 3M 



1/2 



+ 



My 

2r n 3 



-2^ 2 + y 2 



r ° - 3M » + z 2 



r - 2M 



[2M^ 2 (4r G - 15M) 



14r G 5 (l - 2M/r ) 1 /2(r - 3M) 
+ ^ 2 (14r 2 - Q9Mr + 89M 2 ) + 2Z 2 (r G - 2M)(7r G - 2AM)] , 

In terms of the functions defined above, the THZ coordinates (t, x, y, z) are 

x = xcos(QH s ) — ysm(QH s ) 

and 

y = x sm(QH s ) + y cos(QH s ) 
where — Qy^l — 3M/r G , along with 

MZ 



(B5) 

(B6) 
(B7) 



z = rcos(9) + 
+ 



2r 3 (r - 3M) 
MXZ 



[-^ 2 (2r - 3M) + y 2 (r D - 3M) + Z 2 (r a - 2M)] 
[MX 2 (13r - 19M) 



14r G 5 (l - 2M/r ) 1 /2(r - 3M) 
+ ^ 2 (14r 2 - 36r G M + 9M 2 ) + Z 2 (r - 2M)(Ur - 15M)] ) 



(B8) 
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and 

t = t s (l - ZM/r ) 1 ' 2 - ^ — 

V 7 ; (l-2M/r G ) 1 /2 



r 2 (l-2M/r ) 1 /2( ro _3M) 



+ 14r..(r.-2M)(r,,-Mfl " + ^ 

+^ 2 (13r G 2 - 45r Q M + 31M 2 ) + Z 2 (13r G - 5M)(r G - 2M)] . (B9) 

The set of functions (t, x, y, z) forms a non-inertial coordinate system that co-rotates with 
the particle in the sense that the x axis always lines up the center of the black hole and the 
center of the particle, the y axis is always tangent to the spatially circular orbit, and the z 
axis is always orthogonal to the orbital plane. The spatial coordinates are all Lie derived by 
the Killing vector, C^x = C^y = C^z = 0. 

The x a coordinates (t,x,y,z) are locally inertial and non-rotating in the vicinity of T, 
but these same coordinates appear to be rotating when viewed far from T as a consequence 
of Thomas precession as revealed in the QH S dependence in Eqs. (B6) and (B7) above. 

The determination of the THZ coordinates was tedious but not difficult. We looked 
for the relationship between the THZ coordinates and the usual Schwarzschild coordinates 
X A = (t s , r, 9, 0) by using the usual rule for the change in components of a tensor under a 
coordinate transformation, 

9 ~ g W^dx^ (B10) 

where g AB is the Schwarzschild geometry in the Schwarzschild coordinates. The terms 
through 0(X 2 ) (the X in the order term represents any of X, y, or Z) in the definitions 
of t, x, y and z and the rotation represented in Eqs. (B6) and (B7) were chosen so that t 
measured the proper time on the orbit, g a b\r — Vab and d c g a b\r = 0. This much could be 
done easily by hand and resulted in coordinates that satisfied condition (i) in Appendix A. 

The 0(X 3 ) and 0(X A ) terms were found by use of GRTENSOR running under MAPLE 
following a procedure similar to that in Appendix A except that we used homogeneous poly- 
nomials of the form Q 1 i j k X l X 3 X k and ( a ijkiX 1 X 3 X k ' X 1 along with Eq. (BIO) to determine 
C a ij... which resulted in the satisfaction of conditions (ii) and (Hi) in Appendix A. 

Note that ultimately the THZ coordinates (t, x, y, z) are linear combinations of products 
of C°° functions of the Schwarzschild coordinates, and so are C°° functions themselves. 

It is convenient to note that the natural Minkowski metric that goes with the Thorne- 
Hartle coordinates is r] a bdx a dx b = —dt 2 + dx 2 + dy 2 + dz 2 , while its components in the 
original Schwarzschild coordinates are 

dx a dx b , . 

VAB = Vab dxAdxB - (Bll) 

Another form for 7]^ is 

Vab = -V a t V b t + V a x V b x + V a y V b y + V a z V b z (B12) 

or 

Vab = -V a t Vbt + V a x V b x + V a y V b y + V a 5 V b z 

+ (x 2 + y 2 )V a (nh s )Vb(nh s ) + 2[xV [a y - yV (a S:]Vb)(nh s ). (B13) 
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And from the above definitions it readily follows that C^p 2 = C^z = C^V a t = 0, while C^x, 
C^y, and C(t are nonzero. 



APPENDIX C: INTEGRALS OVER $ 

The approach in this and the following Appendix, is similar to that of Appendices C and 
DofRef. [13]. 

In SectionV we define 

X = l-asin 2 $ (CI) 

where 

a = — — — . (C2) 

And we use 

( x -p) = <(l -asin 2 $)- p ) = - f ' (1 -asin 2 $)- p rf$ 

71 Jo 

= 2 F 1 (p,^,l;a) = F p . (C3) 

This result follows almost immediately from 

2 r /2 



2 Z^' 

( X ~ p ) = - / (1 -asin 2 $)- p cM> 

= I [ 1 r 1 / 2 (l-t)- 1 / 2 (l-at)- p dt, 
n Jo 



(C4) 



where t = sin 2 $, and from the integral representation of the hypergeometric function, 
equation (15.3.1) of reference [27] 

2 F 1 (a, b; c; z) = fn^Z b) f ~ W^i 1 ~ tz )~ a dt i ^e(c) > ^e(6) > 0. (C5) 

Two elementary special cases of Eq. (C3) are 

(1 - a sin 2 $) = 2 Fi(-l, ^; 1; a) = 1 - X -a = F_i (C6) 

and 

<(1 - a sin 2 $)^> = 2^(1, 1; 1; a) = (1 - a)" 1 / 2 = F l5 (C7) 

The latter is used in Eq. (63) and leads to the A r term in Eq. (12). The special cases p—\ 
and p — — \ are also easily represented in terms of complete elliptic integrals of the first and 
second kinds respectively, 

-K{a) = 2 F l { l -^l-a) = F l/2 (C8) 
71 2 2 

and 

-E(a) = 2 Fx(-i \, 1; a) = F_ 1/2 . (C9) 

7T Z Z 
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APPENDIX D: LEGENDRE POLYNOMIAL EXPANSIONS 

We require the coefficients A P J 2 {5) in the expansion 

oo 

(5 2 + l-u) p/2 = ^A p e (6)P t (u), for 5^0, (Dl) 
1=0 

for both positive and negative odd-integral values of p. Note that if p is a positive even 
integer then the left hand side is a p/2 degree polynomial in u and the sum terminates with 
I = p. 

First we analyze the negative odd-integral values of p via induction. The generating 
function for Legendre polynomials is 

oo 

(l-2tu + t 2 y 1/2 =^2^P e (u), \t\<l. (D2) 

£=0 

With T defined from 

t = e~ T , (D3) 

Eq. (D2) implies 



(e T + e- T -2 U y 1/2 = J2 e ~ (e+1/2)Tp ^), T>0. (D4) 



-1/2 = 

£=0 

The expansion 



e + e —2 + T+ 0(T 4 ), T -> 0, (D5) 
followed by the substitution 

T = 8V2 (D6) 

in Eq. (D4) provides 

A; 1/2 = V2 + 0(£5), 5^0, (D7) 

which is used in the second term of Eq. (60) to arrive at Eq. (66) and leads to the B r term 
of Eq. (12) and to the absence of a term in Eq. (12) which might have resulted from the 
third term of Eq. (60). 

Differentiation of both sides of Eq. (D4) with respect to T yields 



-\ (e T + e~ T - 2u)- 3/2 (e T - e~ T ) = £ - (t + \) e~^ T P e (u), T > 0. 



(D8) 



Simplification and expansion about T = gives 



(e T + e~ T -2 U y 3/2 = f^^±^P e (u)[l + 0(lT)}, T - 0, (D9) 



2T 



and repeated differentiation extends this result to 



[f + e-r - 2u)- k - 1 ' 2 = ± 2(2 <^ [1 + Q«D] ■ T - 0. (DID) 
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Finally, for k > 1 the expansion and substitution of Eqs. (D5) and (D6) result in 

^ k ~ 1/2 = s^W=T) [1 + om s ^°- (D11) 

For k = 1 

A; 3/2 = ^±1[i + 0(£S)], 5^0, (D12) 

is used in the first term of Eq. (60) to obtain Eq. (63) and, subsequently, the A r term of 
Eq. (12). 

Next, for positive odd-integral values of p in Eq. (Dl), first let p — 1 and multiply the 

left hand side of Eq. (D2) by (1 - u) 1 / 2 and the right hand side by J2e A Y' 2p t( u )- Then > 
integrate over u from —1 to 1; the right hand side is 

J2t £ A 1 / 2 2/(2£+l) (D13) 

e 

from the normalization of the Legendre polynomials, 

1 p t (u)P e (u)du = (D14) 

1/2 

Now, expand the left hand side in powers of t to determine the A/ . This results in 

'A/9 



"-'"'^ i2MP?3) f 'M (D15) 

and 



A '' 2= (2i-i£ + sy < D16) 

The latter is used in the e 1 term of Eq. (60) to obtain Eq. (69) and, subsequently, the D r 
term of Eq. (12). 

For other positive odd-integral values of p > 1 in Eq. (Dl), consider the Legendre poly- 
nomial representation of (1 — w) fc+1 / 2 , with k a positive integer, 

oo 

(1 _ U )*+V2 = J2A k + 1/2 Pe(u) (D17) 

£=0 

which defines the expansion coefficients A^ +1 ^ 2 . The first coefficient Aq +1 ^ 2 is obtained by 
multiplying both sides of Eq. (D17) by 1 = Po(u), integrating over u from —1 to 1 and using 
the orthogonality of the Legendre polynomials to yield 



ofe+l/2 
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The coefficients A\ +l ^ 2 for i > 1 are obtained from Eq. (D15) by induction on k. The 
derivative of Eq. (D17) provides 



-!*+?) f;^ 



-1/2 



fc-1/2; 



£=0 



p/ p/ 

.•-1/2 -^-1 -Q+l 

2£ + 1 



(D19) 



where the prime denotes differentiation with respect to u, and the last equality follows from 
equation (12.23) of reference [28]. A re-summation of this last expression yields 



OO s 1 x OO 



A 



,4 fc-1/2 

2£ + 3 2£ - 1 



fc-1/2 
£+1 



P' 



For £ > 1 



^ +1/2 



k + 



1 



fc-1/2 

e+i 



A 



fc-1/2 



2£ + 3 2^-1 



(D20) 



(D21) 



provides A k t +1 ^ 2 in terms of A\ 1 ^ 2 with the help of Eq. (D18). The final result is 

^fc+1/2 = p k+1/2 (2e+l)/ [(2£ -2k-l)(2£-2k + l)... (2£ + 2k + 1)(2£ + 2k + 3)] , (D22) 
where 

Vk+i/2 = (-l) fc+1 2 fc+3 / 2 [(2A; + l)!!] 2 . (D23) 

for k a positive integer or zero. Eq. (D22) is used to furnish the I dependence in the E\ 
term of Eq. (12). 

The significant conclusions of this Appendix are summarized in Eqs. (D7), (Dll) and 
(D22). 



APPENDIX E: INTEGRATION OF THE SCALAR WAVE EQUATION 

The scalar field resulting from a charge q moving in a circular orbit of the Schwarzschild 
geometry is most easily found following an approach similar to that of Breuer et al. [29] or, 
more recently, Burko [15]. The wave equation for the scalar field is 

V 2 -0 = -Aug, (El) 

where the scalar field source g, being distinct from p, represents a point charge q moving 
through spacetime along a worldline T(r), described by coordinates z a (r). This source is 

g(x) = q J{-g)-^5\x a -z\r))dr 

= q{-g)- l,2 {dt/dT)- l 5\x* - **(*)), (E2) 
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with r the proper time along the worldline. For a circular orbit at radius r Q , expanding g 
in terms of spherical harmonic components provides 



Q = q 



{- 9 y 1/2 5{r - r Q )5{9 - vr/2)5(0 - Qt)S(t - t(r))dr 
r- 2 q5(r - r )5(6 - n/2)5((f) - Qt)dt/dT 

4nr n 



= E f^^ r - ^e^YUO, 0), (E3) 



where 



u m = —mQ, (E4) 
4nq Y; m (n/2,0) 

qtm = JTT3 ' l b5 ) 

r Q at far 



and 



dt 



dr y/l-3M/r ' 

Also, decomposing ip provides 



(E6) 



and the £m component of the scalar wave equation becomes 



djkrn 2(r — M) djj em 
d r 2 r(r — 2M) dr 

We know that 



+ 



uj 2 r 2 £(£+!) 



(r - 2M) 2 r(r - 2M) 



^ = -^^5{r-r ). (E8) 



Yt,- m = (-1)™^, (E9) 
and the reality of g and of the final solution for ip(t, r, 9, 0) requires similar expressions for 
g*_ m and if)t,-m- 

The boundary conditions of interest require only ingoing waves at the event horizon 

il>t m = e™-/r, r^2M, (E10) 
and only outgoing waves at infinity 

il>t m = e- iur '/r r^oo, (Ell) 

where 

r* = r + 2Mlog(r/2M- 1). (E12) 
An expansion of ip£ m starts the numerical integration at large r. We assume that 

*Ur) = e —Y,% (E13) 

n=0 

and, with Eq. (E8), obtain a recursion relation for a„: 

n(n- 1) 1) M(n-l) 2 _ 
an = ^ " -an-i L -an-2, E14 

29 



with the starting values of ao = 1 and a n< o = 0. This is an asymptotic expansion, and we 
begin the integration of Eq. (E8) at a value of r which is just big enough that the sum in 
Eq. (E13) reaches machine accuracy before beginning to diverge. We numerically integrate 
Eq. (E8) in to the radius of the orbit r Q : this provides us with a homogeneous solution ip°° 
with proper boundary conditions at large r. 

Similarly an expansion of ip£ m for small r — 2M starts the numerical integration near the 
event horizon. We assume that 



Y J K{r-2M) n 



(E15) 



n=0 



and, with Eq. (E8), obtain a recursion relation for b n : 

12iwM(n-l) + (2n-3)(n-l) 



(£ 2 + £+V 



12icoM(n 



2M(AinuuM + n 2 ) 
2) + (n-2)(ra-3) 



J n— 1 



£(£+l) i 



AM 2 (Ainu M 
iuj(n — 3) 
2M 2 {AinuM + n 2 ) n ~ 3 ' 



(E16) 



with the starting values of b — 1 and b n<0 = 0. We begin the integration of Eq. (E8) at 
a value of r — 2M which is just small enough that the sum in Eq. (E15) reaches machine 
accuracy within a reasonable number of terms. We numerically integrate Eq. (E8) out to 
the radius of the orbit r Q : this provides us with a homogeneous solution ip H with proper 
boundary conditions near the event horizon. 
The retarded field is 



ret 

£m 



AipP, r<r 



Bip°° 



tm 
oo 
tmi 



r > r n 



with the match at r Q determined by the 5— function source of Eq. (E8), 



B 



dtp 



oo 



dr 



A 



dip, 



P.m. 



dr 



r -2M 



which yields 



and 



A 



B 



,h dip. 



oo 



dr 



1>. 



in i 



oo 
im 



dr 



oo 
tm 



dr 
dr 



= -1>. 



oo 
tm' 



2M' 



Qtr, 



The i component of the radial self-force for ip 



ret 



■IM 

in Eq. (11), is then given by 



(E17) 



(E18) 



(E19) 



(E20) 



-rret 



E 



dip ret 



tm 



dr 



(E21) 



Section VI describes the efficient use of J^f in the determination of the self-force. 



[1] S. Detweiler and B. F. Whiting, Phys. Rev. D 67, 024025 (2003), gr-qc/0202086. 



30 



[2] P. A. M. Dirac, Proc. R. Soc. (London) A167, 148 (1938). 

[3] B. S. DeWitt and R. W. Brehme, Ann. Phys. 9, 220 (1960). 

[4] Y. Mino, M. Sasaki, and T. Tanaka, Phys. Rev. D 55, 3457 (1997). 

[5] T. C. Quinn and R. M. Wald, Phys. Rev. D 56, 3381 (1997). 

[6] T. C. Quinn, Phys. Rev. D 62, 064029 (2000). 

[7] K. S. Thorne and J. B. Hartle, Phys. Rev. D 31, 1815 (1985). 

[8] X.-H. Zhang, Phys. Rev. D 34, 991 (1986). 

[9] Synge, Relativity: The General Theory (North Holland, Amsterdam, 1960). 
[10] L. Barack and A. Ori, Phys. Rev. D 61, 061502(R) (2000). 
[11] L. Barack and A. Ori, Phys. Rev. D 66, 084022 (2002), gr-qc/0204093. 
[12] L. Barack, Y. Mino, H. Nakano, A. Ori, and M. Sasaki, Phys. Rev. Lett. 88, 091101 (2002). 
[13] Y. Mino, H. Nakano, and M. Sasaki, Prog. Theor. Phys. 108, 1039 (2002), gr-qc/0111074. 
[14] L. Barack and A. Ori, Phys. Rev. D 67, 024029 (2003). 
[15] L. M. Burko, Phys. Rev. Lett. 84, 4529 (2000). 
[16] L. Barack and L. M. Burko, Phys. Rev. D 62, 084040 (2000). 

[17] C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravitation (Freeman, San Fransisco, 1973). 
[18] K. S. Thorne and S. J. Kovacs, Astrophys. J. 200, 245 (1975). 

[19] T. Damour and B. R. Iyer, Ann. Inst. H. Poincare (Phys. Theorique) 54, 115 (1991). 

[20] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C: 

The Art of Scientific Computing (Cambridge University Press, Cambridge, 1992), 2nd ed. 
[21] S. Detweiler, Phys. Rev. Lett. 86, 1931 (2001). 
[22] C. O. Lousto, Phys. Rev. Lett. 84, 5251 (2000). 
[23] C. O. Lousto, Class. Quantum Grav. 18, 3989 (2001). 

[24] L. Barack and C. O. Lousto, Phys. Rev. D 66, 061502 (2002), gr-qc/0205043. 
[25] S. Weinberg, Gravitation and Cosmology (Wiley, New York, 1972). 
[26] T. Damour and B. R. Iyer, Phys. Rev. D 43, 3259 (1991). 

[27] M. Abramowitz and I. A. Stegun, eds., Handbook of Mathematical Functions (Dover Publica- 
tions, New York, 1965). 

[28] G. B. Arfken and H. J. Weber, Mathematical Methods for Physicists (Harcourt/Academic 

Press, San Diego, 2001), 5th ed. 
[29] R. A. Breuer, P. L. Chrzanowski, H. G. Hughes III, and C. W. Misner, Phys. Rev. D 8, 4309 

(1973). 

[30] We use ip vet throughout in places where other authors have used ip iuU or ^ total to denote the 
"actual" field [2]. 

[31] The differentiability of ip R is controlled by boundary conditions and initial data. We con- 
sider nondifferentiable initial data or shock waves coming in from boundaries to be physically 
unreasonable. 



31 



